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We study the real time dynamics of the Bose Hubbard model in the presence of time-dependent 
hopping allowing for a finite temperature initial state. We use the Schwinger-Keldysh technique to 
find the real-time strong coupling action for the problem at both zero and finite temperature. This 
action allows for the description of both the superfluid and Mott insulating phases. We use this 
action to obtain dynamical equations for the superfluid order parameter as hopping is tuned in real 
time so that the system crosses the superfluid phase boundary. We find that under a quench in the 
hopping, the system generically enters a metastable state in which the superfluid order parameter 
has an oscillatory time dependence with a finite magnitude, but disappears when averaged over a 
period. We relate our results to recent cold atom experiments. 
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I. INTRODUCTION 

Ultracold atoms trapped in optical lattices (THl] are 
highly versatile systems in which parameters can be 
tuned over wide ranges. The ability to tune these param- 
eters in real time has opened the possibility of studying 
the dynamic traversal of quantum phase transitions ei- 
ther in a "quantum quench" or with a more general time 
dependence. This protocol has received considerable in- 
terest [E|-[l4[ as the resulting systems give examples of out 
of equilibrium dynamics in interacting quantum systems, 
a class of problem that is still not fully understood. 

When bosons are cooled to lie in the lowest Bloch 
band of the periodic potential, their behaviour can be de- 
scribed using the Bose- Hubbard model (BHM) [l5[ . The 
BHM displays a transition between Mott-insulator and 
superfluid phases as the ratio of inter-site hopping J to 
the on-site repulsion U is changed, as has been observed 
experimentally [3, [l6l - [l9| . This transition has been stud- 
ied extensively theoretically and the equilibrium mean 
field solution is well known |2fj|423| . More accurate deter- 
minations u sing quantum Monte Carlo [24l - |28j and series 
expansions [29j verify the qualitative mean field picture 
|3fj| . In addition to cold atoms, there have also been pro- 
posals to realize the BHM in photonic [3l[ and polaritonic 
systems [32j . 

Experimentally there have been investigations of the 
transition from superfluid to Mott insulator or vice versa 
by loading a condensate (or localized atoms) into an opti- 
cal lattice and then increasing or decreasing the depth of 
the optical lattice 0, [U [34| . Both the hopping between 
sites and the on-site interactions in the BHM used to de- 
scribe this situation depend on the strength of the opti- 
cal lattice potential [15| , but the hopping is considerably 
more sensitive to the lattice depth than the interactions. 

Extensive theoretical effort has been expended on try- 
ing to understand the effects of time dependent J/U in 
the BHM (which can allow for a traversal of the phase 



transition). Both sweeps from one phase to another, ei- 
ther gradually or as a quen ch |35t-52| and periodic mod- 
ulations with time (4^,|53 58] similar to experiments in 
Refs. [H, H3| have been considered. A number of pre- 
dictions have been made for these dynamics, including 
the time dependence of the decay of the superfluid order 
parameter for different explicit forms of the time depen- 
dence of J(t) [4(1 |49|; a nd of a wavevector dependent 
timescale for freezing |40l . l43l |49| upon entering the Mott 
phase from the superfluid. Predictions for the transition 
from the Mott phase to superfluid include the genera- 
tion of vortices via the Kibble-Zurek mechanism, and 
scaling of time dependent correlations with the quench 
timescale [4l[. Such scaling (albeit with different expo- 
nents to those predicted in Ref. |41|) was recently ob- 
served in experiments by Chen et al. [33j]. Studies of 
the extended BHM jU and of quenches in the BHM 
[2 IH, Hl| suggest that non-equilibrium states can per- 
sist for considerable times after a quench, especially for 
final states with small values of J/U. In addition to the 
ratio J/U , time dependence of other parameters, such as 
the chemical potential [62j . or even the lattice itself [63[ 
have also been investigated. 

The generation of out-of-equilibrium states from 
sweeps from the superfluid to the insulating phase (or 
vice versa) of the BHM is generic to dynamical traver- 
sal of quantum phase transitions [5T4l4j] and not limited 
to the BHM. Experimentally it is not possible to access 
zero temperature phase transitions, but as the effects of 
such transitions extend to finite temperature, it is in- 
teresting to allow for thermal effects on the quench dy- 
namics. There has been considerable theoretical work 
on the BHM for non-zero temperature [HI, lo4l - l75l ] , but 
most has focused on the equilibrium properties of the 
model - we allow for the effect of temperature in our out- 
of-equilibrium calculation by assuming a thermal initial 
state. 

The approach we take to study the out of equilibrium 
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dynamics of the BHM is to allow J to be a function 
of time with U constant. Our approach is sufficiently 
general to allow for the inclusion of a trapping potential 
and time dependence in parameters other than J. We 
construct a real-time effective action for the BHM using 
a strong coupling approach that can describe physics in 
both the superfhhd and Mott insulating phases. Various 
strong coupling approaches have been proposed to allow 
description of both phases in equilibrium [76M80I | , and we 
generalize the imaginary time approach used in Ref. (76j 
to real time by using the Schwinger-Keldysh formalism. 
Several authors have previously used Schwinger-Keldysh 
or closed time path 81 -[86 1 techniques to study the Bose 
Hubbard model but have not focused on out-of- 

cquilibrium dynamics. 

Given the assumption of time dependent hopping, we 
obtain the effective action within the Schwinger-Keldysh 
formalism. We then obtain the saddle point equations of 
motion, which we are able to simplify to derive a mean 
field equation for the dynamics of the superfluid order 
parameter during a quantum quench from the superfluid 
phase to the insulating phase of the BHM at fixed chem- 
ical potential. We find that generically the solutions we 
obtain correspond to a final metastable state in which the 
superfluid order parameter oscillates with a finite mag- 
nitude, but averages to zero over a period of oscillation. 
We note that the form of the metastable state depends on 
the value of the chemical potential and relate our results 
to work showing that global mass redistribution is im- 
portant for the equilibration of cold atoms in traps after 
a quantum quench [94j |. 

This paper is structured as follows. In Sec. [IT] we derive 
the effective action using the Schwinger-Keldysh/closed 
time path (CTP) technique and in Sec. IIIII we study the 
saddle point equations of motion for order parameter dy- 
namics. In Sec. II VI we conclude and discuss our results. 



II. EFFECTIVE ACTION 



potential. The Hamiltonian 



H = H v - l^N = ^r^2hi(hi - 1) -/z^ftj. 



contains only single site terms, and Hj contains all of 
the hopping terms - we allow for the possibility that the 
hopping amplitude between sites % and j may be time 
dependent. 



A. Schwinger-Keldysh technique 

The Schwinger-Keldysh [8ll . [H| or closed time path 
(CTP) technique [83-86] is an approach that allows a 
description of out of equilibrium or equilibrium quantum 
phenomena within the same formalism. The usual ap- 
proach to finite temperature calculations is to use the 
Matsubara formalism, which is restricted to equilibrium, 
and requires analytic continuation to obtain real time dy- 
namics. The advantage of CTP methods is that the prob- 
lem is formulated in real time so that out of equilibrium 
problems can be tackled and no analytic continuation is 
required - the price to pay is that the number of fields 
in the theory doubles, a second copy of each field propa- 
gates backwards in time. As discussed by e.g. Niemi and 
Semenoff [84| . the notion of time ordering needs to be 
replaced by that of contour ordering in order to calculate 
Green's functions. 
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In this section we discuss the application of the 
Schwinger-Keldysh technique to the Bose Hubbard model 
and derive a strong-coupling effective action for the 
model. The Hamiltonian for the Bose Hubbard model 
takes the form 



FIG. 1: Contour for the Schwinger-Keldysh technique for a 
system with inverse temperature /3. The value of a is arbitrary 
in the interval [0,0] (Ref. [H). 

For a thermal initial state, as we will assume here, the 



generating functional Z factorizes [84|: 



Hbh = — ^ Jij \ a\aj +a]a 



<ij> 
U 



with Ci, C2, C3, and C4 contour segments as illustrated 
in Fig. [TJ The value of < a < (3 is arbitrary 
work with a = for simplicity. 



Hj + H , 



where &i and a] are annihilation and creation operators 
for bosons on site i respectively, n, = a\a,i is the number 
operator, U the interaction strength, and n the chemical 



B. Effective action for the Bose Hubbard model 

We may write a path integral for the generating func- 
tional of the BHM: 
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[Va*][Va]e 



*5bhm [ a ,o] 



(i) 



where a is a bosonic field and we omit source fields and 
set h = 1. The action for the Bose Hubbard model has 
the form 



Sbhm = / dt [a* a (t) (id t ) Tl b a lb {t)] + Sj + Sv, (2) 



where 



dt J2 J H [< a {t)^ b a jb (t) + a* a (t)T% b a ib (t)] 



(3) 



in which the hopping varies as a function of time to cross 
from the superfluid to the Mott Insulating phase. Hence 
we require a formalism that allows for an adequate de- 
scription of both phases. We thus generalize to real time 
the strong coupling method used in imaginary time by 
Sengupta and Dupuis [Zf3]. The advantage of this ap- 
proach, as pointed out in Ref. [76| is that it leads to a 
normalized spectral function, which allows for the calcu- 
lation of the excitation spectrum and momentum distri- 
bution in the superfluid phase, whilst also giving a good 
description of the Mott insulating phase. A similar equi- 
librium effective action based on the Keldysh approach 
was recently obtained in Refs. [90l - l92l |. 

The approach requires two Hubbard-Stratonovich 
transformations. The first of these decouples the hop- 
ping term. We introduce a Hubbard-Stratonovich field ip 
and make use of the identity (derived in Appendix |A"|) 



and Sjj is the action associated with where ai a is the £—*(£* >7+f»7* 

field at site i on contour a, where a = 1 or 2. We use 
notation such that t % is the i th Pauli matrix, acting in 
Keldysh space rather than spin space. ^ G wr jte 

We perform a Keldysh rotation so that 



xe 



01 (t) 

a 2 (t) 



a q (t) \ _ f ( a t (t) 
a c (i) ~ \a 2 (t) 



where a q and a c are the quantum and classical compo- 
nents of the field respectively [92|, I95l-l97j , and 



f 



The effect of this on the action is that r 3 in the 1 , 2 basis 
becomes t 1 in the g, c basis, hence (dropping tildes) 

/OO 
Jij [ a *a{t) T lb a 3b{t) + a* a (t)T^ b a lb (t)] . 



Z = / [Dip*][Vip]e~^ 



with 



jWty'M = / e - l /rftE 1 '0*a(t)^b^b(t)+^a(t)T Q 1 i) Q* i) (t) 



where the average (. . .) is taken with respect to 

/OO 
dtJ2 KaW (id t )T^a ib (t)]+Su. 
-OO 



(5) 



Unlike previous studies of the BHM using closed time 
path techniques [87H93| , we are interested in the problem 



W[ip*,rp] can be used to calculate the 2n point con- 
nected Green's functions G nc for the bosonic field a via: 



Gia 1 ...a„a' 1 ...a' n (^1 ) • • • ) ^1 ) • • • > ^n) 



(-l) n s^w[1p*,tlj) 



tp* = !/' = () 



*(^ 1 )™' r a 1 1 6 1 ■ • • T a„fc„' r a 1 ' 1 b' 1 • • • T a' n b' n (®ibi • • • O lfc „ (tn) ^ (C) • • • 4^ » 

(6) 



where the superscript c indicates a connected function. Note that the connected Green's function vanishes if not all 
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sites are identical. Thus, we may write (similarly to Ref. 76]): 

(—i\ n r°° 



V'fei (*l) • • • ^io„ (*n)V>i< (C) ■ • ■ V'iai 



Xr ai6i • • ■ T a„b„ T a' 1 b' 1 ■ ■ ■ T a' n b' r fi ' i.b l ...b n b' 1 ...b' rl • • • ) *ni • • • ) t n ), 



(7) 



and so 



where 



cn 

°int — 



(-i) T 

(n!) 2 



E 



^a 1 (*l)---^a„(*n)^(4)...^V 1 (*'l) 



XT ai6i ■ ■ • r a„ b„ T a{ b[ ■ • ■ T a' n b' n ^'?^i 1 ...b„b' 1 ...b' n • • • i *n j *li • • • i O- 

Summarizing the effective action to quartic order after the first Hubbard-Stratonovich transformation gives: 
SlsWM = -\J dtY,^ a {t){Jij)- l Tl h i>ib{t) - j dt 1 dt 2 Y,^LA t i)<b 1 G t b 1 b 2 (h,t2)T b 1 2a2 ^ a2 (t2) 

ij ' i 

+ - / dhdhdtsdUj^^aAt^tazfc) T aib 1 T a 2 b2 



(8) 



(9) 



We discuss how the mean field phase boundary at zero 
and finite temperature may be obtained from Eq. (J9j) in 
Appendix El Sengupta and Dupuis [76j observed that 
although the equilibrium action of the form obtained in 
Eq. leads to the correct mean field phase boundary, 
it leads to an unphysical excitation spectrum in the su- 
perfiuid phase. This can be rectified by performing a sec- 
ond Hubbard-Stratonovich transformation (76[. Starting 



from 



(10) 



introduce a field z such that 



so we have 

Z= [Vz*]{Vz]e l f dt ^^ ( - 2Jij ' )z '> a ^ T " bZjb( - t ' ) / [p^*][p^] e i / dt ^J z >*»( t ) r » 1 ^ i6 ( t ' + ^( t ) T ^ 2i6 ( t )]e w ^*''''']. (12) 



As discussed earlier, 

where Sq is the quadratic term 



So = -E/^« 1 ^)< 6l G i6l62 (t l! * 2 R a2 ^ 2 (t 2 ), 



and let 

£ iW(z*,z) _ I [p^*]^j e »SG+i/dtEj<„(t)r a \^b(0+^a(t)^a\^6(*)] e *E~ =2 S i ? lt (V'%V'^ 



/ e */d*E i [<a(«)^a\V'i6(*)+^*„(*)r a 1 () 2 lb (t)]+iE~ = 2'S i » t (V'*,V')\ _ Mty 

\ ' Sg 
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We next perform a cumulant expansion for W(z* , z) and keep only terms in the action that are not "anomalous" (for 
further discussion see Refs. [7y, |98[) to obtain 

Z = [ [Vz*][Vz]e lS '« [z "- z \ 



where in calculating the effective action to quartic order in z, we truncated j^]™ =2 S" nt — > iS? lt , with 
The effective action to quartic order in the z fields is 



oil 



/ dtJ^z* a (t)(2J ij )T^ b z jb (t)+ I dt 1 dt 2 J2^ ai (h)[G ta2ai (t2,h)} l z ia2 {t 2 ) 

ij i 
+ J y dtidt2dt^dt4 U aia2a3a4 (t\, t 2 , *3) ^4)^! (^l) z i*a 2 {^2)zia 3 (^3)^404 (^4), 



(14) 



where 



1 



x |[Ct05ai(*5)*l)] 1 [G ! io6a 2 (*6)*2)] 1 [^io 3 a^ (*3j [^ia 4 ajj (*4) ^)] 

+((a 4 ,t 4 ) <-> (a 3 ,t 3 )) + ((a 6 ,t 6 ) o (a 5 ,t 5 )) + ((a 6) f 6 ) o (a 5l i 5 ); (a 4 ,t 4 ) o (a 3 ,t 3 ))} . 

(15) 

Following the arguments presented in Appendix B of Ref . 98] , it can be shown that the Green's functions for z are 
the same as those for the original field a. 

We note that the following symmetry relations hold for the interaction kernel u from the definition above: 

U a bcd(tl,t2,t3,t4) = Ubacd(h, t 2 , ^3) ti) = Uabdcifl, t 2 , ^3) ti)- 

It can also be seen from the definition in Eq. that 

Gia 5 a 6 a' a' (*5) t&, ^6> = ^ iaaa^aLai (*6> *5> ^6> = ^iasaeaia;, (*5j ^6: ^5: ^g)- 



Similar symmetry relations in the Keldysh structure of 
four point functions were noted in Refs. [9l|, [92|. Hence 
there are only 8 independent components we need to eval- 

uatp- C 2c C 2c C 2c C 2c C 2c C 2c C 2c and 
uouc. ^qqqqj cqqqi u ' qqqd u ' qqco u ccqql ^cqcq' ^qccc £LUU - 

^cccg- The remaining four point function G^ cc = by 
causality Explicit expressions for each of the non- 
trivial components are written down in Appendix [D] We 
will find that for our study of the simplified equations 
of motion away from the degeneracy points of the Mott 
lobes that we will only require G 2 .^ , but the expressions 
we provide in Appendix [D] allow for a more general study 
of the equations of motion than we provide here. 

The mean field phase boundary can be determined 
from the effective action Eq. (fT4)) from the vanishing of 
the coefficient of z*z c by noting that 



(^i&i(*l)$ii> 2 (*2)) = - iT l ia A^ [<3ia 2 ui(*2,tl)] 1 , 



and that the matrix Green's function takes the form 



G(t u t 2 ) = 



o G£(h,t 2 ) 



where Qq , Gq , and Qq are the retarded, Keldysh, and 
advanced Green's functions determined using the single 
site Hamiltonian Hq respectively. These Green's func- 
tions are discussed in more detail in Appendix [Bj We 
can thus obtain 



G^ihM) 



[Go 1 ]* (tiM [QoTM 

[GoY (*1,*2) 
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where 

[g^] R (h,t 2 ) = [g?(h,h)]- x , (i6) 
[a - 1 ] A (*i,* a ) = [gi(ti,t 2 )V\ (17) 

[Q^] K {tiM) = - J dt'dt" [g^tut')]' 1 

xg?(t',t") [g£(t",t 2 )]-\ 

(18) 

which along with g^itx — t 2 ) = £?o (fa — ii), allows one 
obtain the standard equation for the mean field phase 
boundary [Eq. ([BTT]) ]. 

III. EQUATIONS OF MOTION 

We can obtain the equations of motion for the order 
parameter from the saddle point conditions on the action: 



5S c ff _ n SStff _ 

It is helpful to note that 
[GoSiMT 1 = 0; [G.SiMT 1 = [g^ 1 ]* (hM), 

[G qc {ti,fa)]~ l = [g^ifaMY 1 = KifaM)Y\ 

and 

[G^tufa)}- 1 = [g^fa)]- 1 = [$ *(*a,ti)]~\ 

to obtain the equations of motion as 



= 2J ij {t)z jc {t) + I dfa^itM)] z ic (fa)+ dt 2 [g^] (t,fa)z iq (fa) 

J —OG J— OO 

dt 2 dt 3 dt4U qa2 a 3ai {t, t 2 ,t 3 , t 4 )z* a3 (t 2 )z ia3 {h)z iai (u), (19) 



/CO \ f 

dt 2 [g£(t,t 2 )]~ z iq {t 2 ) + - / dfadt 3 dt 4 (t, t 2 ,t 3 , t 4 )z* a2 (fa)z ia3 (t 3 )z ia4 (t 4 ), 



(20) 



with implied summation over a 2 , a 3 and a 4 . The solution of these two equations is rather involved in the general 
case, but the expressions above allow for the description of the spatial and temporal evolution of the superfluid order 
parameter in both the superfluid and Mott insulating phases. By taking appropriate variations of the effective action 
Eq. (|14[) one may also obtain equations of motion for correlations of the z fields. In order to gain some insight into 
the out of equilibrium dynamics of the situation in which the hopping J is time dependent and there is a sweep across 
the boundary of the superfluid, we derive a simplified equation for the dynamics of the superfluid order parameter 
and study its properties numerically below. 



A. Simplified equation of motion 

To investigate the nature of the solutions of the 
equations of motion, we make some simplifications to 
Eqns. (fT9]) and ([20|). We focus on low frequencies and 
long length scales to determine an equation for the mean 
field dynamics of the order parameter. 

We assume that in the limit t — > — oo, the system is in 



the superfluid phase and the hopping J(t) is not changing 
with time. The initial conditions require z\ = z 2 , which 
implies that initially z q = and z c = \f2z\ , where z = z\ 
is the superfluid order parameter. If z q remains small 
under evolution with time then we can focus only the 
equation of motion for z c : Eq. (fT9|) . To see that this is 
indeed the case, we need to note that (see Appendix |B|) 



Go(") = -^^e^^-^Kr + ^J^ + M-C/O+rJ^ + ^-C/^-l))]. 

r=0 

I 

Hence terms involving g^ will only contribute to the low frequency dynamics when fj, ~ Ur for some integer r. 
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These values of \i correspond to the values of chemical po- 
tential where for J = there is degeneracy between Mott 
insulating states with r and r— 1 particles per site. We re- 
strict ourselves to values of the chemical potential away 
from degeneracy, in which case we only need to retain 
terms involving and Q$ . In order for z q to become ap- 
preciable, the term u cccc (t, t 2 , t 3 , t 4 )z* (t 2 )z c (t 3 )z c (t4,) in 
Eq. (|20[) must be appreciable. This term depends on the 
two particle connected Green's function G 2 q c qqq . Similarly 
to Qq , Gggoo om y contributes to low frequency dynam- 



ics when /i ~ Ur for some integer r. We can thus safely 
ignore z q and focus solely on the dynamical equation for 
z c : Eq. (|19p. Taking into account considerations about 
which terms are important for low frequency dynamics as 
we did above, it turns out that for values of the chemical 
potential away from /i ~ Ur, the only connected function 
that we need to evaluate is (^L , which is specified in 
Appendix [Dl Writing z\ = z, we can obtain a simplified 
form of Eq. (|T9l) by first noting that 



dt 2 [g?(t,t 2 )] 1 z{t 2 



2? 



(21) 



which we can expand using 

[Sol _1 («) = K] :l + - £ [OS] 1 + ^ £ K] 1 

leading to 



ui=0 



(22) 



where 



i & 



2 dto 2 ^ 



j=0 



u=0 



Explicit expressions for v, A and k 2 can be easily computed from Eqs. (IB3|) and (IB4[) and are given in Appendix 
[Cl The temperature and chemical potential dependence of these quantities is displayed in Figs. [2] a) - c). The phase 
boundary of the superfluid phase at finite temperature is shown for reference in Fig. [5] d) . We can see that the 
strongest temperature dependence of the parameters is for values of fi/U close to an integer (which we ignore), and 
that both A and k 2 are relatively insensitive to thermal effects over a wide range of \ijU values. The interaction term 
u is most sensitive to temperature and starts to deviate strongly from its zero temperature value by temperatures 
as large as T ~ 0.2L7, which corresponds to the temperature at which there is full melting of the insulating phase 



Using a similar expansion to the one used to derive the 
mean field phase boundary, we can note that after a 
Fourier transform in space 

d 

2Jij(t)z jc (t) -> 2J(t)^cos(fcj-a)z(k,t) 

3=1 

1 



2J(t) 



d--k 2 a 2 
2 



z(k,t), (23) 



for small ka. We will focus on the long wavelength limit 
and ignore terms of order ka. 

We only retain the k = part of the interaction term, 
in keeping with our focus on long wavelength physics, 
and we take the low frequency limit of the interaction 
term by expanding the two particle connected Green's 
function and the retarded and advanced Green's func- 
tions about the u — limit. Recalling from above that 



|z c (t)| 2 = 2|z(t)| 2 , we may approximate the interaction 
term by —u\z\ 2 z, where u is stated in Appendix [Cl and is 
in accord with the static value calculated for equilibrium 
in Ref. [76j]. Thus we have as our approximation to the 
equation of motion: 



\2dJ(t) + v\ z(t) - iX^ 1 - k 2 ^—^- - u\z{t)\ 2 z{t) = 0. 
at ot z 

Take Jit) = J + j(t), where Jo is chosen so that 
2dJ + v = 0, 

i.e. Jo is chosen to lie on the mean field phase boundary 
for the superfluid for a given fi. Hence we may write the 
approximate mean field equation of motion as 
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FIG. 2: Plots of a) A; b) k 2 ; and c) u as a function of fi/U and inverse temperature /3. d) The phase boundary of the superfiuid 
state is shown in as a function of fi/U and 2dJ/U at several temperatures for reference. The filling per site in the Mott 
insulating phase at zero temperature is indicated. 



r2 "^2" + + S ( t ) Z + U \ Z \ 2z = > 



Of 



(24) 



where S(t) = —2dj(t). Even after the simplifications 
made above, this equation for the dynamics of the or- 
der parameter is a non-linear second order differential 
equation, for which we are not able to find analytic solu- 
tions in general. Below we discuss numerical solutions of 
this equation, along with an analytic solution that can be 
determined in a special case which illuminates the prop- 
erties of the solutions of the equation. 

We study Eq. (|2"4")l for fixed \i and time- varying J. In 
experiment, there is a confining potential so that there is 
a position dependent local chemical potential 

/"local (r) = A* - V(r), 

where V{r) is the trapping potential. The solutions we 
obtain for the dynamics at fixed /x should be compared 
to the experimental situation in which one views the dy- 
namics at fixed radius in a symmetric trap. (This pic- 
ture should be reasonable at time scales shorter than the 



timescale for global mass redistribution in the trap, which 
can be quite long compared to microscopic timescales 

If we fix [a, then there are two possibilities for the dy- 
namics that we should consider: a) the particle-hole sym- 
metric case, in which case A = 0, and b) the generic case, 
in which A ^ 0. The particle-hole symmetric case corre- 
sponds to the transition at the tip of the Mott lobe as 
illustrated in Fig. [3] 

We consider traversal of the quantum critical region as 
S(t) varies with t. We demand that 



lim S(t) 

t— > — oo 



-Sq; and 



lim 8(t) 

t— ¥QG 



Si. 



In our numerical solutions we use the form 

Si - 8 

JQ, 



S(t) 



So + Si \ . .ft 
tann 



(25) 



where, similarly to Cucchietti et al. [4l|, who studied the 
transition from Mott insulator to superfiuid in the one 
dimensional BHM, we assume that there is a timescale 
tq which is the characteristic time for S(t) to cross from 
— <5q to Si. 
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Mott 
Insulator 



Superfluid 



Particle hole 
symmetric case 
X = 



Generic case 



FIG. 3: Two possibilities for quantum phase transition at 
constant /i. 



B. Particle-hole symmetric case 

In the particle-hole symmetric case, A = and the 
saddle point equation takes the form 



n 2 — v + S(t)z + u\z\ 2 z = 0. 



(26) 



We can choose z(t) = p(i)e lS '*\ and then real and imag- 
inary parts of the equation give 

= p-p(o) + 5p + up 3 , 

= 26p + 6p 

where we rescaled 

t = Kt; 5{t) = 5{k£) = S(t), 

and wrote the equations in terms of the rescaled time 
co-ordinate i. The second equation can be integrated to 
give 



i.e. Bp 2 = c, so 



ln^ = -21np + ci, 



p - -7T + Sp + up 3 = 0. 



The initial condition that the system is deep in the super- 
fluid phase implies that as t —¥ -co, — > 0, and p — > 0, 

so p = \J~^i and c = (we choose 9 — without loss of 
generality). Thus 

p + Sp + up 3 = 0. 
Rescaling p — > p/y/u, then dropping the tilde and bar, 
p + 6(t)p + p 3 = 0, 



with p —> V^o as t — >• — oo. In the long time limit, when 
S(t) = St, then we may rewrite the differential equation 
for p as 



d_ 

di 



\iP? 



Si 2 



1 



0. 



Then 



6 lP 2 



and, writing p = £y, t — r/x, we have 



dy 



-g) = {l-k 2 )-(l-2k 2 )y 2 -k 2 y\ 



with 



k 2 - -f 2 n 2 



1 - 2k 2 = Sxr) 2 



and we can solve to get 



k = 



V2 



which must satisfy < k < 1. The solution to our equa- 
tion as t — > oo is thus 



P = £cn 



.V2k' 

which in the original variables is 

Z(t) = ^ Cn (^^ fc 

In general we cannot determine the value of £ analyti- 
cally. We can obtain an analytical solution if there is 
a jump in S(t) from — Sq to +5\ at t = 0. [Note that 
this form of S(t) violates the assumption that we made 
in deriving the equation that frequencies are low, but 
the solution in this case is still instructive, as it shares 
many features with the solution for more physical forms 

of 5(t).] We know z(t) — \J~^- for t < 0, and recalling 

cn(0; k) = 1, we get £ = y/~3~o, which implies 

1 1 



V2 



1 + T 



and so we get 



z(t) 



1 



/ VW+h)t 1 
— cn I ;— - 

U \ K V2 1+ 5i 



which is periodic in time with average value and period 
1 1 



AK 



V2 + / V^o + Si 
V So / 
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FIG. 4: Dynamics of z(t) normalized to unity in the particle 
hole symmetric case, for a variety of tq. The parameters are 
PU = 100, fj, = 0.4142136, k 2 = 0.707107, u = 0.1038, and 
we take So — 1.83, Si — 0.17 = Jo(/x). This corresponds to a 
quench from 2dJ/U = 2.0 to 2dJ/U = 0.0. The inset shows 
the value of Zmax(rQ) as a function of tq. 
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FIG. 5: Dynamics of \z(t)\ normalized to unity in the generic 
case, for a variety of tq. The parameters are j3U = 100, 
p = 0.25, A = -0.28, n 2 = 1.55, u = 0.1277, and we take 
So = 1.85, Si — 0.15 = Jo(/i). This corresponds to a quench 
from 2dJ/U = 2.0 to 2dJ/U = 0.0. 



where 



K(k) 



dO 



\/l - k 2 sin 2 



We obtain numerical solutions of Eq. (|26p with 8(t) 
taking the form given in Eq. (j2"5")l at the particle-hole 
symmetric point in the first Mott lobe for several different 
values of tq, as displayed in Fig. |U One can see that in 
each case, for large values of t tq, the form of the 
solution is that z(t) oscillates in a periodic manner with 
a magnitude that decreases with increasing tq. When 
averaged over a period T at times t 3> tq, 



1 

T 



t+T 



dtz{t) = 0, 



as we would expect in the Mott insulating state. Dehn- 
mg z max (t q ) = lim^oo \z(t)\ we can see that z m ^{T Q ) 
decreases with increasing tq without any indication of 
saturation, as illustrated in the inset to Fig. 2J Note that 
in our numerical simulations t is measured in units of 



C. Generic case 

In the generic case in which A ^ 0, we start with 
Eq. (|24|) and try for a solution of the form 



z(t) =p(t)e ieit \ 

Taking real and imaginary parts of the equation, after 
substitution gives 

p - p (V) 2 \ - X P 6 + 5(t)p + up 3 = 0, 

k 2 (29p + Op) + Xp = 0. 



Integrating the second equation with respect to t leads 
to 



A „2 



2 2 
K Z p Z 



(27) 



In the t — > — oo limit, p and are constant, so we can 

A<5 



determine c = -|/9q 
equation for p: 



K 2 p - 



4k 2 p 3 



, , and we obtain the following 



— + S(t)p + up d 



0. 



(28) 



We solve Eq. {28} numerically for a variety of values of 
tq and display \z(t)\ = p(t) for p/U = 0.25 (well away 
from both degeneracy and the particle hole symmetric 
case) in Fig. [5j 

The solution displays the similar feature to the 
particle-hole symmetric case that the average of z over a 
period (z) T = 0. However, it is clear that as tq increases, 
there does not seem to be any decay in the values of |z(t)|. 
By rescaling the time with tq , we can see that in fact the 
different traces collapse onto each other, as we display in 
Fig. El 



D. Chemical potential and temperature 
dependence of dynamics 

The traces of z(t) and \z(t)\ that we displayed in 
Figs. HE] were for a particular value of the chemical po- 
tential in the generic case and for a low temperature 
((3U = 100) in both cases. It is of interest to see whether 
the observation that in the non particle-hole symmetric 
case that there is a metastable state after a quantum 
quench is robust to variations of chemical potential and 



temperature. Defining 



lim 
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FIG. 6: Dynamics of normalized to unity in the generic 
case, for a variety of tq with time rescaled by tq. The pa- 
rameters are as in Fig. [5] 



calculated this for 0.1 < fi/U < 0.9 and temperatures 
ranging from f3U = 100 to j3U — 2. We focus only on 
the first Mott lobe, but from perusal of the chemical po- 
tential and temperature dependence of the parameters 
A, k 2 and u in Fig. [2] we expect that similar qualitative 
results should be obtained for other Mott lobes. We find 
that apart from the particle- hole symmetric point, where 
we believe the displayed finite value of z msK is an arte- 
fact of our numerical calculations, that the transition to 
a metastable state in which z max ^ is generic for a wide 
range of values of fi and persists to temperatures compa- 
rable to the melting temperature of the insulator as illus- 
trated in Fig. [7] It should be noted that the physics that 
we have left out of dynamical equation, namely spatial 
dependence of z and also higher frequency components 
of z will presumably lead to equilibration of z at long 
enough times, but as we argue in Sec. IIV) it may well 
be reasonable to expect that the behaviour we identify 
at the mean field level to be experimentally relevant on 
appropriate timescales. 



IV. DISCUSSION AND CONCLUSIONS 

In this paper we have derived a real time effective ac- 
tion for the Bose Hubbard model using the Schwinger- 
Keldysh technique, generalizing previous work that ob- 
tained an equilibrium effective action [76| . This action 
allows for a description of the properties of both the su- 
perfluid and Mott insulating phases. Hence we are able 
to study the out of equilibrium dynamics as the param- 
eters in the Hamiltonian are changed so that the ground 
state is tuned from one phase to another. We obtain the 
saddle point equations of motion and by focusing on low 
frequency, long wavelength dynamics are able to obtain 
an equation of motion for the superfluid order parame- 
ter. We have focused on this case as the simplest example 
of dynamics, but we emphasise that our approach leads 
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in the large tq limit at several temper- 



to equations of motion that can be used to study high 
frequencies and spatial variations of the order parameter 
and its correlations. 

We study the equations of motion by varying the hop- 
ping parameter J as a function of time at fixed chemical 
potential to sweep from deep in the superfluid phase to 
deep in the Mott insulating phase over a timescale of 
order tq. We study the tq dependence of the super- 
fluid order parameter numerically and find that in the 
long tq limit the system generically reaches a state in 
which the time averaged value of the order parameter is 
zero (as would be expected in equilibrium for a Mott in- 
sulator), but the absolute value of the order parameter 
is non-zero. The magnitude of the order parameter in 
the long tq limit appears to vanish only at the particle- 
hole symmetric value of the chemical potential, and grows 
with distance from the particle-hole symmetric value of 
li. The generic final state is clearly an out-of-equilibrium 
metastable state, with equilibration only possibly for the 
particle-hole symmetric case. The generically non-zero 
value of \z(t)\ in the final state indicates that the system 
retains memory of the initial superfluid state, a feature 
which is observed in quantum revival experiments (99l — 
Il0l| that indicate quantum coherence remains even after 
a quench into the insulating phase. 

There have been several other recent theoretical works 
on the out-of-equilibrium dynamics of the Bose Hub- 
bard model that see evidence of the system entering a 
metastable state after a sweep from the superfluid phase 
to the Mott insulating state. Schiitzhold el al. [!(| stud- 
ied the dynamics in the limit of large number of bosons 
per site and found a slow decay of the superfluid fraction 
for a slow sweep from the superfluid phase to the Mott 
insulating phase. Kollath et al. [HJ investigated the one 
and two dimensional BHM numerically with the number 
of bosons fixed to an average of one boson per site and 
found that for small enough values of the final value of the 
hopping, the system reached a non-thermal steady state 
which was relatively insensitive to the details of the initial 
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state. These authors determined whether the system was 
thermal or not by investigating real-space correlations, so 
it is not possible to make a direct comparison with our 
results here. Most recently Sciolla and Biroli [l3| con- 
sidered the infinite dimensional Bose Hubbard model at 
integer filling and also found that the final state after a 
quantum quench of U showed a non-zero superfluid order 
parameter. Similar features have also been reported for 
mean field studies of fermions after a quantum quench 
102]. 

Whilst the emergence of a metastable state after a 
quench from the superfluid to the insulating state is also 
seen in our work, we study a different situation to the 
previous works. We consider a spatially uniform BHM, 
as do Refs. [lU |4(| EH, but we consider fixed chemical 
potential rather than fixed particle number. To compare 
theoretical descriptions of the out of equilibrium dynam- 
ics of the Bose Hubbard model and experiments on the 
quench dynamics of a fixed number of cold atoms in an 
optical lattice, the physical meaning of working with fixed 
chemical potential needs to be discussed. The presence of 
a spatially non-uniform trapping potential means that in- 
stead of viewing the system as having a uniform chemical 
potential, it is often more convenient to view the system 
as having a spatially dependent local chemical potential: 
Miocai(r) = fi—V(r), where V(r) is the trapping potential. 
For a symmetric trap, this implies that a reasonable de- 
scription of the phase the system is in at radius r can be 
determined by using /j,\ oca i(r) - this implies the "wedding 
cake" structure seen in many experiments. Our study of 
the equations of motion at fixed chemical potential would 
then correspond to studying the dynamics of atoms in a 
trap at fixed radius (albeit with radii corresponding to 
certain values of the chemical potential excluded due to 
the approximations we made in deriving the equation of 
motion). 

This viewpoint appears to be borne out in recen t ex- 
periments [33l. l34l[l03l . ll04| and theoretical work (941 1 105 
on quantum quenches for cold bosons. Natu et al. [94 
argue that the very large differences in re laxa tion times 
observed in Refs. [34[ (of order ms) and [l03j (of order 
~ Is) can be understood if one looks at mass transport 
during equilibration. If the average number of particles 
per site remains the same in crossing from the superfluid 
to an insulator, then equilibration can be quick as in 
Ref. [34[ , but if the average number of particles per site 
needs to change, then there must be mass transport and 



the equilibration is slow as in Ref. [103J. The results we 
find for the long time limit of 2 max illustrated in Fig. [7] 
are in accord with this idea. For the chemical potential 
associated with particle hole symmetry, the value of z max 
decays to (close to) zero, whereas for other values of fi, 
Zm&x can be an appreciable fraction of the value of \z\ 
in the initial state. At the particle hole symmetric fi, 
the average number of bosons per site does not change in 
crossing from the superfluid to the Mott insulator [2fJ , in 
accord with the condition for local equilibration without 
mass transport [94j . Global mass transport is not cap- 
tured within our simplified equation of motion, and there 
is no decay of the metastable state and equilibration on 
a longer time scales. 

The main results of our work and their connection to 
existing experimental and theoretical work in the field of 
cold atoms are outlined above, but there are a number of 
future directions that it might be interesting to pursue 
based on what we have done here. First, a more thor- 
ough study of the solutions of the equations of motion 
allowing for spatial fluctuations and higher frequencies 
than we consider here might lead to further insight into 
the dynamics of the Bose Hubbard model. The inclusion 
of a trapping potential would also allow for additional 
contact with experiment [106| . Second, it would be in- 
teresting to add the effects of dissipation [EH, Il07| | , which 
has been shown to renormalize the phase boundaries in 
the BHM. For cold atoms the effects of dissipation can 
probably be ignored, but in ot her realizations of the BHM 
this may not be feasible 1081 ] - 



Recent experimental advances which allow for high 
spat ial resolution in cold atom experiments [U, Il04 109 
112] suggests that there will be advanced capabilities for 
probing the out of equilibrium dynamics spatially as well 
as temporally, suggesting that there are exciting times 
ahead for studies of out of equilibrium dynamics of Bose 
Hubbard systems. 
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Appendix A: Hubbard-Stratonovich transformation 



Starting from the identities (where z = x + iy) 
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it is easy to show that 

f I" dxdy i| z |2 +i ( z * a+za *) ia * a f j" dxdy _i| z | 2 +i( z * a + za *) 

■ ' ' (-in) 

Using these results we may write (with z\ = X\ + iyi, 22 — x-i + iyi and 23 = £3 + iy 3 ) 

dx ^ [ dx ^ [ ^^ e i|^i 2 -i|^i 2 -^3i 2 e i(, 1 (r+^)+^(£+r ? )+^r+^e+ Z 3^+. 3 ^) 

in J (-in) J (-in) 

d Xl d yi dj 2 rfj/2 rfg 3 dj/3 c «i^ 1 | 2 -«|£ 2 -g 1 | 2 -i|g 3 -g 1 | a c »fg;g+5 a f +s^»j+^ 3TJ *-) 
in {—in) (—in) 

where we change variables to Z\ = Zi, £2 = z\ + z%, and 23 = 21+23. After integrating out z\, we get 

<ev+in*) = [ dx 2 d yz dx 3 dy 3 c 2i(i 2 x 3 +y 2 v 3 ) c i(z*$+z 2 e +z 3V - +z; v ) 
(—in) in 

©(z 2 ,22p(^3,2|)e i(z 2 Z3+Z2Z 3) e i (^S+^e , '+^'7-+^) ) ( A1 ) 



where 



dxdy — dxdy 
V(z,z) = — — ; V(z,z ) = . 

in (— in) 



Appendix B: Mean field phase boundary 

One way to determine the mean field phase boundary 
between the superfluid and Mott insulating phases is to 
determine when the coefficient of the quadratic term in 
the action Eq. (j9|) vanishes. In order to do this it is 
helpful to note that 



These expressions can be evaluated using the interaction 
representation 

a\t') — e i (Hu-IJ.fi)t '-t e -i(Hv-[ih)t' 
d(t) = e i(.Hu-»ri)t &e -i(Hu-»ft)t 



we obtain 



(Bl) 

where 0* and Q* are the retarded advanced and _ {Er _^ t i(Er _ 1 _ Mr _ 1))(4 _ i(Er _, r)t > 

Keldysh propagators respectively, with the subscript '^o \ L i L ) — ' e e e 

indicating that these are the propagators associated with 
Hq. The definitions of the propagators are: 



iQ$(t-t') = ig<(t,t') + ig>(t,t'), 
iQ^(t-f) = 6(t-t')[ig>(t,t')-ig<(t,t')} 
iGo(t-t') = e(f -t)[ig<(t,t')-ig>(t,t')} 



with 



ig$(t,t') 



Tr{tf(t)a(t)po} 
Z 

Tr{a(f)at(t')p } 



where we recalled a\r) — \/r\r — 1) and 
\/r + 1 |r + 1). At temperature T, 



iS < (M') 



r=0 '' 



!J i(/ J -i7(r-l))(i-t') f ,-/3(Br-A'r) 



Y°° e -P(E r -nr) 



(B2) 



J 



Hence we have that the retarded Green's function takes the form 
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G R {t u t 2 ) = -i6(t 1 -t 2 ) 



[( r + l) e »(f*-^r)(ti-t 2 ) _ re »(M-t/(r-l))(ti-t2) 



e — } , (B3) 



which simplifies at T — to 



Go{ti,t 2 ) = -<ff(ti - t 2 ) \(n + l) e ^- Un ^^-^ - noe i(^-c/(no-i))(ti-t 2 ) 
For future reference it will also be convenient to note that 



-P{E r -pr) 



r=0 



which simplifies at T = to 

G^(h,t 2 ) = -i [(n + l) e *(f-^no)(ti-t2) 



(r + l) e i (M-fr)(t 1 -t 2 ) + r . e i( M -C7(r-l))(ti-*2) 



(B4) 
(B5) 



+ n e 



»(A»-t/(n -l))(ti-t 2 ) 



Recalling that we can treat this as a single site problem we have 

p = e -/»[¥WA-i)-MA)]. Z = Tr {p } = ^ e -^-Mr) ) 

and -E r = ^r(r — 1). n = N/M, where N is the number of bosons and M the number of sitesa u is determined 
implicitly from 

re -P(Er-v.r) 



n = 



y-°° e -0(E r -^r) 



At T = 0, the value of /x/t/ sets the occupation number, no which takes an integer value r for r — 1 < < r, 

with degeneracies at [i/U = 0,1,2, 

When we Fourier transform the quadratic part of S^ s in space and time we get: 

-/^£E^ :(w ' k)T ^ (w ' k) "/ ^EC^k^.G^K^J^k). (B6) 
We choose the hopping amplitude J i2 - to take the form 



T — / ^ +j(*)' *, J nearest neighbours 
■^W-j 0, otherwise 



for which (with a the lattice spacing) 



Mt) = [J +j(t)]^ cos (M 

'd-h 2 a 2 ) [Jo+j(t)], 



assuming that h « 1. 

Setting j(t) = for now, when we take the cj, k — > limit we can locate the phase boundary by noting when the 
coefficient of the ipqi^c term in the action vanishes: 

! + g^iuj = o) = o. 



2dJ 



Note that the retarded propagator 



oo 

Z-<r=0 



r=0 



(r + 1) 



u-£/r + w + i0 (j,-U(r-l)+u + i0 
I 



(B7) 



at finite T and for T = 

,R, x "0 + 1 



CM = 



n 



u — [7n + u; + i0 /U — U(n — 1) + uj + iO ' 

(B8) 



The advanced propagator may be obtained from 
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and at T = the Keldysh propagator is 

Qq(u>) = -2in[(n + l)S(u + (i- Un ) 
+n 5(cu + a - U(n — 1))] 



(B9) 



At zero temperature we obtain the standard mean field 
equation for the phase boundary between the Mott insu- 
lator and superfluid phases: 



1 



(no + 1) 



no 



2d Jo [i — Uno fx—U (no — 1) 
This may also be expressed as 



= 0. 



(2n + l)-J±Jl- J(2n + 1) + ,P 



(BIO) 



J 



for no > 1 and //+ = — J if no = 0, where J = 2dJ/U and 
|S = jtx/f/. This affirms that the effective action correctly 
predicts the mean field phase boundary at zero temper- 
ature. At finite temperature the corresponding equation 
is 



1 



2dJn Z 



1 oo 



r=0 



1 



[i-Ur n-U(r-\) 



= 0. 



(Bll) 



The phase boundary as determined from this equation for 
a variety of j3 values is displayed in Fig. H]d). This phase 
boundary is the edge of the superfluid phase - the Mott 
insulator is strictly defined only at T = 0, and at non- 
zero temperature there can be a normal phase separating 
superfluid and insulator, with full melting of the insulator 
for T* ~ 0.2 U [H, [zO]. 



Appendix C: Parameters in the equation of motion 



There are three parameters that enter the equation of motion: 



= K]:io-^=-^K] 



. K 2 1 9 2 r R] -l 

= 2d^ [G ^ 



These can be evaluated to give 



L^t—O 



(r+1) r 

fi-Ur jU-C/(r-l) 



(CI) 



z ^ 

r=0 



(r + 1) r 
[fi-Ur] 2 ~ [/i - U(r - l)] 2 



(C2) 



o A V ■r-^ 
K = > I 

v Z ^ 



(r + 1) 



[fx - Ur} 3 [fi - U(r - l)f 



(C3) 



and 



4 oo 



2Z 



r=0 



4(p + l)(p + 2) 



4p(p- 1) 



4(p+l) : 

[C/p- A i] 2 [2 A1 -(2p+l)C7] T p7(p-l)- M ]2[fJ(2p-3)-2 M ] [/i-ETp] 

4p 2 4p(p + 1) 4p(p + 1) 

~ [t/(p - 1) - M ] 3 ~ [tf(p - 1) - n] 2 \M - t/rf ~ [U( P - 1) - Mn^ - E/p] 2 



The expressions for v, A and k 2 simplify somewhat in the zero temperature limit: 



(y-Un )(p-U(no-l)) 



A = 



(2710-1)^-2^1 , (ii-Uno)(ji-U(no-l)) 



(v + uy 



and 



K 2 = 



2n U-fi f (Ufx(2n + 1) - U 2 (2n 2 - 1) 



(M + C/) 2 V 



(a* + ur 



(C4) 



(C5) 
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Appendix D: Evaluation of the four point function 



To evaluate the four time correlation functions, there are several basic correlations we need: 



a (ti, 12,13,14) 



I Tr {e-^^-"*)a(ti)a(t2)at(t 3 )at(t4)} (Dl) 

OO 

^ — ' 

1 00 

+ l) 2 e^' Bp_Aip ^* 1+ * 3 ~ t2_ * 4+i ' 9)+ ^- Ep+1_Al(p+1 ^ (t2+ * 4 ~* 1_ * 3 ) (D2) 

Z — ' 

00 

— ^p(p + l) e i ( B p-w)(ti+t3-t2-i4+i/9)+i(Bp-i-M(p-l))(t4-t 3 )+i(-E p+1 - AJ (p+l))(t 2 -t 1 )^ pg-j 

— 00 

— y^P(P+ l)e i(Bp ~ w)( * 1+ * 3 ~* 2 ~* 4+i ^ +i(Bp+1 ~^ (p+ ^ (D4) 

p=0 

1 00 

— y^p 2 e '(- E p-W)(*l+«3-t2-t4+tf)+»(-Ep-l-^(p-l))(t2+t4-tl-t3) i 

^ p=0 

_ 00 

S atataa (ti,f 2 ,i3,i4) = — 2)e^ Bp_w ^ tl_ * 4+i ^^ + ^ Bp " 1_M ^ _1 ^^ 2+ * 4_tl_t3 ^ + ^- Ep " 2_A1 ^ _2 ^^ 3_ * 2 ^ ) . (D6) 



tS (ti, t 2 , 13,I4j 



a [ti, t 2 , t 3 , 14J 



i> (ti, t 2 , 13, C4J 



nd^aa^/i + + + \ 

a [t\, 12,13,14) 



p=0 



In addition we require the two point correlations 



1 00 

= _ y^(p + i) e »(S p - W )(ti-t2+i/3)+*(£p+i-M(P+l))(*2-ti) 



W(ti,t 2 ), (D7) 



C° to (ti,t 2 ) = 

p= 

*^(t 2 ,ti)- (D8) 



pe i(Bp- w )(ti-t 2 +i/3)+i(Bp-i- ([1 (p-l))(t2-ti) 



The actual expressions are rather tiresome to derive but are given here for completeness, where we use the notation 
9ij = 9(ti —tj): 
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Gqqqq (*1>*2, Hi ti) 



{ [#12#23 + 021014] (V W (i 3 , t 2 , t!,t 4 ) + B° W (t 4 , *1,* 2 , t 3 )) 



2 

2 

+ [012024 + 021013 
+ [013032 + 031014 
+ [013034 + 031012 
+ [014042 + 041013 
+ [014043 + 041012 
+ [032021 + 023034 
+ [042021 + 024043 
+ [023031 + 032024 
+ [043031 + 034042 
+ [024041 + 042023 



5 aW {tsMMM) + B oW (*4, i 2 ,ti,t 3 ) 

ryaa^ aa^ /< , , , \ , rya^aa^a/, , , , \ 
O (t 2 , 13) £l,£4 J + -D (14, 11,13) 12 J 

ryaaa^a^ (• 1 1 * \ , Tya" aaa" ( t 1 1 1 \ 
£> (12, tl, £3, £4 J + -D ( > £4,£3 J £l J £2j 

ryaa^aa*' (+ + + + \ \ r>a^aa^a(-L + + + \ 
D (t 2 , £4, £l, £3j + -D (£3, £l, £4, £2j 

ryaaa^a^ /> > > > \ , rya* a* aa ( ± > > > \ 
tS (t2, £1, £4, £3j + -D (,£3 ) £4,£l,£2 / ) 

jyaaa^ 1 • 1 1 1 \ , rya^a^aa/j. ± ± ± \ 

E> (,£i, £2, £3, £4; + -d (,£4,£3j£2,£iJ 

ryaaa^a^ /> > > > \ , na^ a 1 aa ( ± > > > \ 

a {11,12,14,13) + a \H,t4:,t2,H) 

ryaa^aa^ /> > > > \ , rya* aa* a ( > > > > \ 
O \tl,ts,t2,t4) + D \t4,t2,t 3 ,ti) 

Tyaa)a*a[+ + + + \ \ jyaa* a* a f 1 + + + \ 

a (,£i, £3, £4, £2j + -d (£2,£4,£3,£iJ 



B aa aa (£1, £4, £2, £3) + B a aa a (£ 3 , £2, £4, £1) 

ryaa^ a / 1 > > > \ , ryaa* a* a ( ± > > > \ 

d (ti, £4, £3, £2j + ts (£2,£3,£4,£iJ 



+ [034041 + 043032 

-i { [C aaf (h,t 3 ) + C aia (t 3 , £x)] [C aat (£ 2 , £4) + C ata (£ 4 , £ 2 )] 
+ [C aaf (£ 1; £4) + C ata (£ 4 , £1)] [C aaf (£ 2 , £3) + C ata (£ 3 , £ 2 )] } , 



(D9) 



Glg qq {tl,t 2 ,h,t4) 



B aaa a (£i,£ 2 ,£ 4 ,£ 3 ) - B a a aa (£ 3 ,£ 4 ,£ 2 ,£i) 

ryaa^aa^ (1 » 1 1 \ Tya" aa" a ( t 1 1 1 \ 
-D (,£l; £3j £2, £4j — -D (,£4, £2, £37 £lj 

ryaa^a^a/i » 1 < \ jyaa^a^a/i > > > \ 

ti (ti, 13,14,12) — a (12,14,13,11) 

ryaa^aa^ /+ < > > \ rya* aa* a ( ± > > > \ 
ti (tl,l4,t2,t3) — tf [t3, t 2 , £4, £lj 

J^aa) O* a l ■ ■ • • ^ ^nn^n^n 



2 | — 021 [032 + 023034] ^B aaa a (£1, £2, £3, £4) — B a a aa (£4, £3^2, £l)^ 
-021 [042 + 024043] 

— 031 [023 + 032024] 
-031 [043 + 034042] 

— 041 [024 + 042023] 

— 041 [034 + 043032] \B aa a °(£i, £4, £3, £2) — B aa a a (£ 2 , £3, £4, £l) 

— [021013 — 031012024] (yB a aaa (£3, £l, £2, £4) — B a aaa (£4, £ 2 , £1, £3) 
+ [021014 — 041012023] aaa (£3, £2, £l, £4) — -B" aaa (£4, £l, £2, £3) 

— [031012 — 021013034] (^B aaa a (£ 2 , £1, £3, £4) ~ B a a ""(£4, £3, £1, £2) 
+ [031014 — 041013032] B aa ' aa ' (£ 2 , £3, £1, £4) - B a ' aa ' a (£ 4 , £1, £3, £2) 

— [041012 — 021014043] (^B aaa a (£2, £l, £4, £3) — B a a aa (t 3 , £4, £1, £ 2 ) 
+ [041013 — 031014042] ^£> aa aa (£2, £4, £l, £3) — -B" aa °(£3i £l, £4, £2)) j 

-i {031 [C ata (£ 3 , £1) - C aa " (ti, £3)] [C aat (£ 2 , £ 4 ) + C a+a (£ 4 , £ 2 )] 



+041 



ata (£ 4 , £1) - C aaf {h, £ 4 )] [C aat (£ 2 , £3) + C ata (£ 3 , £ 2 )] } 



(D10) 
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G 2 c c C qq{ti-,t2,h,ti) — - ^0 32 9 2 i (^B aaa a (ti, t 2 , t 3l ti) + B a a aa (i 4 , t 3 , t 2l tl)j 

+ 042021 

+S 3 i[Si 2 9 23 — 8 32 9 2i ] (^B aa aa (t±, t 3 , t 2 , ti) + B a aa a {t4,t 2 ,t 3l ti)j 

+041 [032024 — 042023 

)] (i?- W ( tl ,U, t 2 ,t 3 ) + B aWa (i 3 , t 2 , U, h)) 

— 031042 (B ootot °(ti, t 3) U, t 2 ) + B aa ' a ' a (t 2 ,U, i 3 , ti)) 

— 041032 (V aVa (t 1; u, t 3 ,t 2 ) + B aa ' a ' a (t 2 ,t 3 , u, hfj 

+ 031012 

^oaatat + £T Vaa (t 4 , t 3 , t U t 2 )) 

+0 4 i0i2 ^S aaa a (i 2 ,ii,i 4 , £3) + B a a aa {t 3 ,ti, t\, t 2 )j 

+ 032 [041013 — 031#14] 

(V aW (t 2) t 3) t u u) + B aWa (i 4 , hMM)) 

+ 042 [031014 — 041013] 

+ [031012024 + 042021013] ^£> aaaa (f 3 , ^ , t 2 , + B a aaa {ti,t 2 ,t\,t 3 )) 
+ [041012023 + 032021014] (^B a (t 3 ,t 2 , t\,ti) + -B a a<m (f4, t\, t 2 , t 3 )) j 

-2 {031042 (C aa+ (t 2) t 4 ) - C a ' a (t 4 , t 2 )) (C a ' a (t 3 , h) - C aat (il,t 3 )) 

+ 041032 (C aa+ (t 2 , t 3 ) - C ata (t 3 , i 2 )) (C ata (t 4 , *l) - C aat (t l9 1 4 )) } , (Dll) 



G 2 c c qcq {ti,t 2 ,t 3 ,ti) 



^{021 [043032 — 023034] (^B aaa a (t\ , t 2 , t 3 , t 4 ) + B a a aa {ti,t 3 ,t 2 ,tif) 
+ 041 [023034 — 043032] ^£> aaa a (t\,ti,t 3 ,t 2 ) + B aa a (f 2 , f 3 , t 4 , 1 1 )^ 
+ 043 [021014 — 041012] ^B aaa (t 2 , t\, ti, t 3 ) + B a ° aa (t 3 ,t i ,ti,t 2 )) 
+ 023 [041012 — 021014] 

£ B aW (t 3 ,t2,h,U) + B a W (t 4 , *1, t 2 , t 3 )) 

5 aaaV {t u t 2 , U, t 3 ) + B ata " aa (t 3 ,t 4 , t 2 , h) 



— 021043 
+ 023031 
+ 043031 

— 023041 
+ 041013 
+ 021013 

+ [021013034 + 043031012] (^B aaa a {t 2l ti,t 3l ti) + B a a aa {ti , t 3 , t\ , t 2 )j 
+ [023031014 + 041013032] (fl" oW {t 2 ,t 3 , h,U) + B a ' aa ' ° (t 4 , h , t 3 , t 2 )) } 
-Z023041 fc ata (i 4 ,tl)-C aQt (t 1 ,i 4 )l \c aa \t 2 ,t 3 )-C a ' a {t 3 ,t 2 ) \ , 



nflfl T aa T / . . . . \ . na^aa^a/i . • • \ 

£> (11,13,12,14,) + a (,14,12,13,11; 

i> (11,13514,12; + -D (12,14, 13, ti; 

a (11,14,12,13; + -0 (13, t2, 14, ti; 
a (12,145 11,13; + -D (13,11,14,12; 
B ataaat {UMMM) + S ataaat (t 4 , t 2 , ti, is) 



(D12) 
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Gll cq (ti,t 2 ,t 3 ,U) 



2 |^43^32^21 (yB a a aa (ti,t 3 ,t 2 ,tl) 
+043042021 (B aaa " a \t U t 2 MM)- 

— 042023031 (B aa ' aa \ti,t 3 ,t 2 ,U)- 

+043031042 

— 041042023 ^-B a<1 aa (tl, *4, ^) — 

— 043032041 ( y B aa a a (ti, ti, t 3 , t 2 ) — 

— 043031012 (^B aaa a (t 2 ,ti,t 3 ,ti) — 
+041012043 

+041013032 

— 041013042 ^-B aa aa (^2,^4,^1,^3) — 

— 042021013 aaa (i 3 , ti, i 2 , ^4) — 

— 041012023 (V W (i 3 ,*2,ii,t 4 )- 



B aaatat (t l7 t 2 ,t 3 ,t 4 )) 

(^3; ^4, ^2, ^l) 



^ a a'aa 



B 



a ' aa ' a 



(^4 j ^2, ^3, *l) 



naa T o T o j. j. j. \ 
Dd^aa^a/-t + + + \ 

ti it3,£2,Z4,ZiJ 

jDaa^a^a(4. + + 4. \ 
13 \t2,t3,t4,ti) 

r»a^ a* aa ( + + + + \ 
ti (£4jC3jC1jC2; 

Tja^a^aafs. ± ± ± \ 

n (13^4,11,12) 

Tia) aa) a ( + , , , \ 
ti \t4,ti,t3,t<2,) 

TDaJaaJa/j. + + + \ 
ti (13,11,14,12) 

jja^ aaa* (± + j. + \ 
ti (t4,t2)tl,l3) 

Sa^ aaa^ (+ j. j. + \ 
(I4, %i,%2, 13) 



}■ 



(D13) 



Gqqqdti, t 2 , t 3 , ti) ~ - |^34 [023 + #32#2l] (^B aaa ° (t\ , t 2 , t 3 , ^4 ) — B" ° 00 (^4 , ^3 , h , i 1 

+ [#24043 - ^34042021] (V aaV [tl,t 2 ,U, h) - B a ' a ' aa {t 3 , ti, t 2 ,tl)) 
+024 [032 + 02303l] (^B aa aa (ti,t 3 , t 2 , £4) — B" ° (*4 , t 2 , *3 , h 
+ [034042 - 024043031] (fl 00 * * ^!, t 3 , ti, t 2 ) - B aa ' a ' a (t 2 , ti, t 3 , tl)) 
+ [014042023 — 02404l] ^-B"" °° (tl, t4, *2, £3) — B a aa ° (t 3 , t 2 , ti, tl )) 
+ [014043032 — 03404l] (^B aa a a (ii , 1 4 , t 3 , t 2 ) — B aa a a (t 2 , t 3 , ti,i )) 
+034 [013 + 031012] (^B aaa a (t 2 , ti,t 3 , ti) — B a a aa (ti, t 3 , ti, t 2 )) 
+ [014043 — 034041012] (^B aaa a (t 2 ,ti,ti,t 3 ) — B a a aa (t 3 , ti, t\, t 2 )j 
+014 [031 + 013032] (^B aa aa (t 2 , t 3 , ti,ti) — B a aa a (t4,ti,t 3 ,t 2 )j 
+ [024041013 — 014042] ( y B <ia °° (t 2 ,ti, tl,t 3 ) — B a aa a (t 3 , t\,ti, t 2 )j 
+024 [012 + 021013] 

(> W (i 3 , ti, t 2 , ti) - B aW (ti, t 2 , h, t 3 )) 
+014 [021 + 012023] (yB a aaa (t 3 , t 2 , tl, ti) — B a aaa (ti,tl,t 2 ,t 3 )^ j 

-i {024 (C aat (ti, t 3 ) + C a ' a (t 3 , h)) (c aat (t 2) t 4 ) - C ata (t 4 , t 2 )) 
+0i 4 (c aat (i 2 ,i 3 ) + C ata (i 3 ,t 2 )) (C aat (i 1 ,t 4 )-C ata (t 4 ,ii))}, (D14) 
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G 2 q c qcc {ti,t 2 ,h,U) 



\ {# 23 #34 (B aaa ' a ' (ti,t 2j t 3 , U) + B atataa (t 4 , UMM)) 

+#24#43 

+#24 [#13#32 ~ #23#3l] (-B° a (*1 , *3 > ^2, £4) + B a aa a {ti,t 2 ,t 3 ,t 1 ) 
+ [#13#34#42 + #24#43#3l] (^B aa a a (t\, t 3 , ti, t 2 ) + B aa a a {t2,t4,t 3 ,ti)j 
+#23 [#14#42 — #24#4l] ^-B aa aa (tl, *4, ^3) + B a °° "(£3, <2, ti, tlfj 
+ [#14#43#32 + #23#34#4l] ^-B aa (il , £4, ^3 , ^2) + B aa a a (t2, t 3 , ti, tifj 
+#13#34 

^aaatat ^ ^ ^ ^ + B ototoo (t 4 , t 3 , *1, * 2 )) 
+014043 (^B aaa a (t 2l t\, ti, t 3 ) + B a a aa (t3,ti,ti,t 2 )j 
+#14 [#23#31 — #13#32] 
+#13 [#24#41 — #14#42] 

^aaW ^ ^ ^ ^ + B „t oa t a ^ ^ ^ ^ 

— #13#24 

— #14#23 (-6° aaa (£3, *2, t\, ti) + _B a aaa (ti,ti,t2,t 3 ) S j^ 



+« #13#24 



C a a (t 3 , h) - C aal (h,t 3 ) C a ' a (ti, t 2 ) - c aal (t 2 , u) 



G q C ccc {tl,t2,t 3 ,ti) 



+#14#23 [C a " a (ti, h) - C aat (t U ti)] [C a " a (t 3 , t 2 ) - C aa " (t 2 , t 3 )] } , 
2 j#12#23#34 ^B aaa a (ti,t 2 ,t 3 ,ti) — B a a aa (ti,t 3 ,t2,ti)^j 



(D15) 



+ #12#24#43 


naaa^a^ ( + + + + \ 

a \t\,t2,ti,t 3 ) 


jja^a^aaii + + + \\ 
— tS [t 3 , ti, 12, t\ ) j 




+#13#32#24 


iz>aa!* aa* ( > > < > \ 
O \t\,t 3 ,l2,ti) 


-na*aa^a/< > > < \\ 

— a \pi,t2,t 3 ,t\) \ 




+ #13#34#42 


T3aa*a*a(t > < < \ 
& \t\,t 3 ,ti,l2) 


— B aa a a (t2,ti,t 3 ,ti)j 




+#14#42#23 


jz>aa* aa* ( + + + + \ 
D \t\,ti,l2,t 3 ) 


r>a^aa^a( + + + + \\ 

— a \t 3 ,t2,ti,ti) j 




+#14#43#32 


Tjaa!* a!* a ( > > > > \ 
D \t\,ti,t 3 ,l2) 


Tjaa^a^a/i > > > \ \ 
— tS \t2,t 3 , 14,11 J 1 




— #12#13#34 


Tjaaa) a) ( > > < < \ 
& (t2,ti,t 3 ,ti) 


-na*a^aa{< . . . \\ 
— ti \pi,t 3 ,t\,l2) \ 




— #12#14#43 


nftflfl* a** ft 1 . . \ 

a \t2,t\,ti,t 3 ) 


T^a*a)aas + + + + \\ 
— ti (£3, t4,ti, 12 ) j 




— #14#13#32 


j^aa) aa* 1 > > > > \ 
a (t2,t 3 ,ti,ti) 


— ti \pi,t\,t 3 ,l2)\ 




+#13#14#42 


Tjaa" aa" (> > < < \ 

a \t2,ti,t\,t 3 ) 


r>a'aa'a{, > > > \ \ 
— ti [t 3 , ti,t\, 12 ) J 




— #13#12#24 


jz>a* aaa* ( + + + + \ 
& \t 3 ,t\,l2,ti) 


Tja^aaa** (> > > + \\ 
— ti \pi,t2,t\,t 3 ) j 




+#14#12#23 


( iz>a* aaa* ( > > > < \ 
Itf (t 3 ,t2,ti,tij 


-B afaaa \ti,t 1 ,t2,t 3 ))Y 


(D16) 
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